Van der Waals loops and the melting transition in two dimensions 
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Evidence for the existence of van der Waals loops in pressure p versus volume v plots has for 
some time supported the belief that melting in two dimensions is a first order phase transition. 
We report rather accurate equilibrium p(v) curves for systems of hard disks obtained from long 
Monte Carlo simulations. These curves, obtained in the constant volume ensemble, using periodic 
boundary conditions, exhibit well defined van der Waals loops. We illustrate their existence for 
finite systems that are known to undergo a continuous transition in the thermodynamic limit. To 
this end, we obtain magnetization m versus applied field curves from Monte Carlo simulations of 
the 2D Ising model, in the constant m ensemble, at the critical point. Whether van der Waals loops 
for disk systems behave in the L — > oo limit as they do for the 2D Ising model at the critical point 
cannot be ruled out. Thus, the often made claim that melting in 2D is a first order phase transition, 
based on the evidence that van der Waals loops exist, is not sound. 
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Unphysical looking "loops" in pressure versus volume 
curves have been coming out of approximate calculations 
for nearly a century [EL These so called van der Waals 
loops have also been showing up in computer simulations 
of melting for over three decades |^-0[. As Mayer and 
Wood pointed out M, pressures that increase with vol- 
ume, that would be ruled out by van Hove's theorem for 
macroscopic systems Q , are indeed to be expected when 
simulating melting of finite systems, van der Waals loops 
that decrease as system sizes increase have been observed 
in simulations ^|J7) . Their existence has almost invariably 
been taken as evidence of a first order phase transition 
[^r(|,[§ (though not always ||) and this has contributed 
much to the often held belief that the solid-fluid phase 
transition in two dimensions (2D) is first order pp| , pT| . 

The purpose of this paper is threefold: (1) to give ex- 
amples of van der Waals loops that do sometimes show 
up for finite systems that undergo continuous phase tran- 
sitions in the thermodynamic limit; (2) to point out that 
since their size (defined below) is exactly equal to the 
free energy barrier for nucleation of the other phase, it 
follows that van der Waals loops are to be taken as signs 
of first order transitions only if their size vanishes in the 
thermodinamic limit as the inverse of the linear system 
size (L); (3) to report accurate data for van der Waals 
loops that we have obtained for two dimensional systems 
of 256 and 1024 classical hard disks, in the fixed volume 
ensemble, and to show that their size dependence is in 
very good agreement with more extensive data that fol- 
low from simulations in the constant pressure ensemble 
(also known as the NpT ensemble) [12 13] that seem to 



point to a second order transition, rather than to a first 
order one. 

The pressure p(v) exerted by a system with a given 
fixed volume per particle v is usually obtained from 
Monte Carlo (MC) or from Molecular Dynamics simu- 
lations carried out at constant volume. In order to ob- 
tain p(v) one makes use of expressions that are derived 
from the virial theorem fl4|| , which in turn follow from 
the relation 



p{v) 



df(v,T) 
dv 



(0.1) 



where T is the temperature and / is the Helmholtz free 
energy per particle. We have performed long Monte 
Carlo simulations (1.2 x 10 s MC sweeps in each run, of 
which the first 0.3 x 10 s sweeps are allowed for equilibra- 
tion) in the canonical ensemble for systems of 256 and 
1024 hard disks. The results obtained are shown in Fig 
1 (as o and □ for N — 256 and 1024, respectively). 

Throughout this paper, "volume" v actually stands for 
the area of a two dimensional system; it is given in terms 
of the closest packing area vq, and has therefore no units. 
The pressure p is actually a force per unit length, which 
we give in terms of kT/vQ and has therefore no units. 

There is an alternative way to obtain the same func- 
tion, p(v), that illustrates how van der Waals loops come 
about for finite systems. Consider the probability den- 
sity, P p (v), that a system at a given pressure p have spe- 
cific volume v. P p (v) can be obtained through Monte 
Carlo simulations carried out at a given pressure p, in 
the NpT ensemble. Data for P p (v) that have been ob- 
tained |lij for a system of 256 hard disks in the solid 



1 



and fluid phases are shown in Fig. 2a. Data for P p (v) 
exhibiting coexistence of both phases are shown in Fig. 
2b [Q. The p(v) curve that ensues in the canonical en- 
semble (that is, —df/dv), can be obtained from P p (v), 
since P p (v) and f(v) are related by 



P p (v) oc exp{-[/0) +pv]N/kT}, 



(0.2) 



where k is Boltzmann's constant. Each data point (v,p) 
exhibited in Fig. 1 cither as a • or as a ■, and in Fig. 3 
as a •, has been obtained from an independent MC run at 
the corresponding value of p. Note that any two volumes, 
such as v a and Vb in Fig. 3, that fulfill p(v a ) = p{vb) are 
most probable volumes (see Fig. 2) when p, instead of v, 
is fixed. On the other hand, v c , which is the portion of the 
loop where dp/dv > and satisfies p{v c ) — p{v a ) = p(vb), 
is the least probable volume. 

Alternatively, f(v) may, of course, be extracted from 
(at least in principle) P p (v) obtained from a single simu- 
lation at an arbitrary constant p, using Eq. (2). Pressure 
versus volume curves follow from Eq. (1). Data points 
thus obtained, from P p (v) for p = 7.64 and N = 256 (ex- 
hibited in Fig. 2b), are shown in Fig. 3 (as □). The good 
agreement between independent sets of the data points 
in Figs. 1 and 3 gives an indication of the accuracy of 
our equilibrium results. 

For comparison, we also plot in Fig. 3 p versus the 
mean volume (v) that is obtained in the NpT ensem- 
ble for a system of 256 disks (shown as O). No van der 
Waals loops obtain. This is because, in the NpT en- 
semble, d(v)/dp = —N((v — (v)) 2 )/kT, which is clearly 
negative. 

We next give an example that underscores the fact 
that while van der Waals loops follow for finite systems 
from first order phase transitions, the converse is not true. 
Consider the Ising model in 2D. The applied field h versus 
magnetization, to, is shown in Fig. 4 for Ising systems of 
64, 256 and 1024 spins, for periodic boundary conditions, 
at the critical temperature T = T C . These curves are un- 
usual because we obtained them from MC simulations in 
an constant m ensemble, rather than in the more often 
used ensemble in which h is fixed. In our simulations, 
we keep to constant using the Kawasaki algorithm [ jEs) . 
In it, one spin is flipped up while another spin is flipped 
down at each MC step. The canonical and NpT ensem- 
bles, discussed above, correspond to the constant to and 
constant h ensembles, respectively. 

We arrive at the h(m) values for the constant magneti- 
zation ensemble given in Fig. 4 as follows. Consider first 
a simulation performed at constant h, in which spins are 
flipped, and the total magnetization M (given by Nm) 
is consequently not conserved. Let w(m + 2/N <— m) be 
the conditional probability that a system known to have 
an M value make the transition M + 2 <— M when no ex- 
ternal field is applied. Since Ph{m) oc exp[— Nf(m)/kT] 
for h = 0, we can write the detailed balance condition, 



w(m + 2/N <— m) 



2/N) 



-Nf(m+2/N)/kT 
e -Nf{m)/kT • 



Let the Hamiltonian TL become H — hM when field h 
is applied. Now, h plays no role in a calculation in the 
constant M ensemble, but an h(m) can be obtained for 
the given value m from the relation h = df(m)/dm. 
Taking logarithms of both sides of the above equation 
gives f(m + 2/N) — f(m) in terms of the transition rates. 
The approximation f (m + 2 /N) - f(m) = (2 /N) [df(m + 
1) /dm] gives h{m + l) = (N/2)[f(m + 2/N) -f{m)]. Wc 
thus arrive at 



kT 

h ( m ) = — ln 



w{m + l/N <- to - 1 /N) 
w(m-l/N + 1 /N) 



(0.4) 



after shifting to — > to — 1 for symmetry's sake. In or- 
der to obtain the transition rates, we proceed as follows. 
First note that the probability for an up spin flip from 
a given spin configuration is proportional to either 1 or 
exp(— AE/kT), depending on whether the corresponding 
energy change AE is either negative or positive, respec- 
tively. Accordingly, after each MC sweep, having applied 
Kawasaki's rule throughout the entire system, we assign 
to each spin down either the number 1 or the number 
exp(— AE/kT), if flipping it up would lower its energy or 
raise it by AE, respectively. (No spin is actually flipped.) 
The sum of such numbers [1 and exp(— AE/kT)] over all 
down spins in the system averaged over an MC run is our 
unnormalized estimate of w(m + l/N <— m — l/N). 

Alternatively, the same curves for h(m) may be ob- 
tained from simulations in which h is fixed, by applying 
the relation h = df(m)/dm to probability curves Ph(m) 
that follow from such simulations. Results obtained in 
this fashion, from both the constant m and the constant 
h ensembles, are exhibited in Fig. 4. 

We have also obtained h(m) curves (not shown) for 
T < T c for the 2D Ising model. However, loop sizes 
for T < T c and for T = T c vary rather differently with 
system size. By loop size, we mean the free energy 
Ag = J h(m)dm over the domain defined by to < and 
h > 0. Data points for LAg/kT are shown in Fig. 5 for 
T/T c = 0.9, 0.95 and 1. Data points for the Gibbs free 
energy 



\p{v a ) -p(v)]dv 



(0.5) 



(0.3) 



(the shaded area in Fig. ||) are also shown in Fig. [| for 
disk systems of various sizes. In order to make Ag unique, 
we choose v a such that p(v a ) is the Maxwell construction 
pressure, p m Q, that is, j^\p(v) - p m ]dv = 0. ^From 
data points shown in Figs. 1 and 3, we obtain, making 
use of Eq. (3), the two data points shown in Fig. 5 as • 
for systems of L x L spins for L = 16 and L — 32. Data 
from previous simulations in the NpT ensemble are also 
shown (as o) for L = 16, 20, 24 and 32 (l3). The error 
bars shown in Fig. 5 follow from the procedure described 
in the Appendix. 

As Lee and Kosterlitz have explained in some detail 
Q , the macroscopic limit L Ag (a measure of the surface 



2 



tension) does not vanish for first order phase transitions. 
This follows from t he fo llowing simple argument. As may 
be seen from Eq. fl0.5| ) and comparison of Figs. || and 
||, L d Ag (the system's dimension d is 2 in this case) is 
the free energy barrier that is surmounted by the system 
when a fluctuation takes it over the top, at v c , starting 
from the (locally) most probable volume v a - We may 
therefore think of L d Ag as the free energy of the wall 
that arises between two coexisting phases when one of 
them is nucleated from the other one in order to make 
a transition from one phase to the other one. Since the 
wall thickness is finite for first order transitions, it follows 
then that L d Ag ~ which is the desired result. 

Inspection of Fig. 5 shows that vanishing of LAg in 
the L — ► oo limit, as for the 2D Ising model at the crit- 
ical point, can clearly not be ruled out for disk systems. 
Thus, the often made claim that melting in 2D is a first 
order phase transition, based on the evidence that van 
der Waals loops exist is not sound. Further re- 

sults for larger systems would help to establish how the 
L — ► oo limit of the surface tension behaves. 

It is perhaps worth stating explicitly that whereas 
phase coexistence and nonvanishing surface tension that 
are associated with first order transitions imply van der 
Waals loops, their appearance depends on boundary con- 
ditions when continuous phase transitions are involved. 
Indeed, whereas Ph(m) for the 2D Ising model at the crit- 
ical point is bimodal for periodic boundary conditions, it 
exhibits one single maximum (no van der Waals loops 
then) for free boundary conditions p|. Analogously, no 
van der Waals loops are obtained for systems of disks for 
hard crystalline walls Jl3| for N < 4096 or for systems 
of disks on spherical surfaces JT9j ] . This provides support 
for the proposition that melting of systems of hard disks 
in 2D unfolds through a continuous transition. 
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APPENDIX A: 

We specify here how we obtain the error bars shown in 
Figs. 1 and 5. 

The error bars shown in Fig. 1 for data that follows 
from MC runs in the canonical ensemble are obtained as 
follows. We divide each MC run into 5 "time" intervals, 
and calculate an average pressure value for each one of 
the five intervals. Twice the values of the standard devi- 
ations thus obtained from each such set are exhibited in 
Fig. 1 as error bars. 

There are two kinds of error bars shown in Fig. 5. 
We first discuss the ones that follow from MC runs in 
the NpT ensemble. We have obtained the probabilities 



P p {v) from independent MC runs for various values of p. 
As discussed in the text, the free energy f(v) is given by, 



-Nf(v)/kT 



P P (v] 



Npv/kT 



(Al) 



Using six curves for f(v)/kT thus obtained, we exhibit 
N[f(v) +p m v]/kT in Fig. 6, for various values of p and 
p m = 7.865, for systems of 1024 disks. Slightly different 
values of LAg/kT are obtained from each one of those 
f(v) + p m v curves (see Fig. 6), from which the standard 
deviation is obtained. It is shown as the error for L = 32 
in Fig. 5. Other error bars shown in Fig. 5 for the val- 
ues of LAg/kT that follow from the NpT ensemble are 
obtained similarly. 

We next explain how we obtain the error bars for the 
data points for LAg/kT shown as • in Fig. 5. These 
errors follow from the errors shown in Fig. 1 for values of 
p(v) that were obtained from simulations in the canon- 
ical ensemble. An explanation is called for because the 
error evaluation procedure is not trivial: variations in 
p(v) lead to variations in the Maxwell construction pres- 
sure and in the integration limits of Eq. (3.5). However, 
a bit of reflection shows that such changes in the lim- 
its of integration give contributions to errors in Ag that 
are of second order in Sp(v) [where Sp(v) is the differ- 
ence between an erroneous pressure and the correct one] . 
Accordingly, we shall neglect variations in v a , v c , and 
Vb (which are defined as in Fig. 1). Thus, the first order 
change, Sp m , that is induced in the Maxwell construction 
pressure p m , by errors in p(v), is given by, 



(Pm + Sp m )(v b - V a ) 



\p(v) + Sp{v)}dv, (A2) 



which, by virtue of the definition of p m itself, leads to, 



Sp-n 



(Vb -V a ) J v 



5p(v)dv. 



(A3) 



Now, the first order variation S Ag, which follows from 
Eq. (0.5), is given by, 



S(Ag) = 5p m (v c - v a ) 



Sp(v)dv, 



(A4) 



which, upon substitution of p m from Eq. (A3), becomes 

2 r v b rv c 

$(Ag) = t r[(v c - v a ) / Sp(v)dv - (v b - v c ) / 6p(v)dv] 

( v b ~ V a ) J Vc J Va 

(A5) 

Finally, we replace integrals by sums. Furthermore, we 
note that all errors obtained for p(v) for different values 
of v follow from different MC runs and are therefore as- 
sumed to be statistically independent. We thus arrive at 
the average value of [S(Ag)] 2 , 



(IHAg)} 2 ) 



2\l/2 _ 



Av 



(Vb - V a ) 



(q + s), 



(A6) 



3 



where 



q=(v c -v a )[ £ ^) 2 ] 1/2 , 



s=(v b -v c )[ £ 5pK) 2 ] 1/2 }, 
and Av — Vi+i — Vi. This is the desired expression, 



(A7) 
(A8) 
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FIG. 1. Pressure p versus volume v data points from MC simulations systems of N disks. The "volume" v stands for an area, 
it is given in terms of the closest packing area vo. The pressure p is a force per unit length, given in terms of kT/vo. Neither v 
nor p have therefore any units, o and □ stand for results from simulations in the constant volume ensemble for N = 256 and 
N = 1024, respectively. • and ■ stand for results, extracted from simulations in the NpT ensemble making use of Eqs. (1) 
and (2), for N = 256 and iV = 1024, respectively. All data points follow from runs of approximately 10 8 MC sweeps, after 
equilibrating the system for 3 x 10 7 MC sweeps. For details about the error bars shown, see the Appendix. 
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FIG. 2. (a) Frecuency of occurrence P P (v) for specific volume v for a system of N = 256 for p — 7.83 (solid line) and for 
p = 7.40 (dashed line). The units for v and p are given in the caption for Fig. 1. (b) Same as for (a) but for p = 7.64. These 
curves follow from runs of over 2 x 10 8 MC sweeps. Lines shown go through datapoints obtained, one for each Av = 10 -3 bin. 
Volume values where d[f(v) +pv]/dv — are marked with arrows. 



G 




FIG. 3. Data points for p vs volume for systems of N = 256 disks. O stand for data points (< v >,p) obtained from 
simulations using the NpT ensemble. The units for v and p are given in the caption for Fig. 1. □ are for the numerically 
obtained derivative —df(v)/dv from the frecuency of occurrence P p (v) for a single pressure, p = 7.64, value. • are for points 
(po,v) fulfilling the relation d[f(v) +pov]/dv = where f(v) +pov follows from ln[P p (v)] for p = po- For example, for p = 7.64 
(marked with a dashed line in the figure) we find three different solutions (v a ,v c ,Vb) (marked with arrows as in fig. 2b). 




FIG. 4. Magnetic field h versus magnetization m for the 2D Ising systems of L x L spins at the critical temperature, for 
L — 8, 16 and 32. Continuous lines stand for data that follows from probability Ph(m) curves, obtained from simulations in 
the constant h ensemble for h = 0. The umbrella method was used to obtain P h (m), covering the whole range of — 1 < m < 1 
values with 16 "umbrellas". Each one of the three continuous lines shown follows from 16 MC runs of 10 8 sweeps over the 
entire system for N > 12 and 5 x 10 6 sweeps for L — 8. ■, o, and O stand for L = 8, 16 and 32 respectively, obtained in the 
conserved magnetization m ensemble (with the Kawasaki algorithm). Every point shown follows from a MC run of 5 x 10 6 
sweeps for L = 8, and 10 s sweeps for N > 12. 
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FIG. 5. Data points for LAg/kT for a system of L x L disks. • and O stand for data obtained from systems of L x L hard 
disks in the constant volume and NpT ensembles, respectively. O, □ and + stand for the 2D Ising model at temperatures 
T/T c = 0.9,0.95, 1, respectively. Lines are only guides to the eye. Error bars for data points shown as O are approximately 
given by the size of the circles, except that the error is (±0.0006) approximately 4 times smaller for L — 32. For a detailed 
account about how these errors, as well as the error bars shown for the • data points, were obtained, see the Appendix. 
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FIG. 6. Quantity N[f(v) + p m v]/kT versus v, up to a constant, obtained from MC simulations of systems of 1024 disks for 
each of the values of p shown; p m = 7.865. Each one of the six curves shown follows from a MC run of at least 2 x 10 8 MC 
sweeps over the whole system. Ideally, all curves should be equal. We define the free energy barrier, AG as the local maxima 
of N[f(v) + p m v] minus the average value of its two minimum values. All values of AG/kT thus obtained are shown in the 
figure. We obtain the average value AG/kT — 0.62. The corresponding standard deviation is 0.02. 
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